#include <stdio.h>
#include <stdlib.h>
#include <math.h> 
#include "system.h"
#include "ran_uniform.h"
 
// generates initial positions/velocities
void Init(void)
{
  int i,j,k,Number,Nplace;
  double tmp,Size;
  double Uold=0,EMStep,ForceValue; 
  
  
  // generate velocities from a Gaussian and calculate the kinetic energy 
  UKinetic=0.0;
  for(i=0;i<NumberOfParticles;i++)
  {
    Velocities[i].x=RandomGaussianNumber();
    Velocities[i].y=RandomGaussianNumber();
    Velocities[i].z=RandomGaussianNumber();
    UKinetic+= 0.5*SQR(Velocities[i].x)+SQR(Velocities[i].y)+SQR(Velocities[i].z);
  }
  
  tmp=sqrt(Temperature*3.0*NumberOfParticles/(2*UKinetic));
  for(i=0;i<NumberOfParticles;i++)  {
    Velocities[i].x*=tmp;
    Velocities[i].y*=tmp;
    Velocities[i].z*=tmp;
  }
  UKinetic *=tmp*tmp;
 
  // calculate initial positions on a lattice
  Number=(int)(pow(NumberOfParticles,1.0/3.0)+1.5);
  Size=Box/(Number+2.0);
  Nplace=0;
  for(i=0;i<Number;i++) {
     for(j=0;j<Number;j++) {
       for(k=0;k<Number;k++)     {
         if(Nplace<NumberOfParticles)   {
           Positions[Nplace].x=(i+0.01*(RandomNumber()-0.5))*Size;
           Positions[Nplace].y=(j+0.01*(RandomNumber()-0.5))*Size;
           Positions[Nplace].z=(k+0.01*(RandomNumber()-0.5))*Size;
         }
         Nplace++;
       }
     }
  }
  
  
 //Energy Minimization: steepest descent algorithm
 EMStep = 0.01*Box;
 for(i=0; i<500; i++)  {
	if(i==0)  {
		Force();
		Uold = UPotential;
		printf("Initial Potential:%10.4f\n",Uold);
	}
	
	for (j=0; j<NumberOfParticles; j++) {
		ForceValue = sqrt(SQR(Forces[j].x)+SQR(Forces[j].y)+SQR(Forces[j].z));
		Positions[j].x += EMStep*Forces[j].x/ForceValue;
		Positions[j].y += EMStep*Forces[j].y/ForceValue;
		Positions[j].z += EMStep*Forces[j].z/ForceValue;
		//periodic boundary condition
		if(Positions[j].x<0)  
			Positions[j].x += Box;
		else if(Positions[j].x>=Box)
			Positions[j].x -= Box;
		if(Positions[j].y<0)  
			Positions[j].y += Box;
		else if(Positions[j].y>=Box)
			Positions[j].y -= Box;
		if(Positions[j].z<0)  
			Positions[j].z += Box;
		else if(Positions[j].z>=Box)
			Positions[j].z -= Box;
	}
		
    // calculate new potential energy
    Force();
    
    // check if the new positions are acceptable
    if (fabs((UPotential-Uold)/Uold) < 0.0001) {
	    printf("Equilibriumed at step: %4d, UKinetic/N: %8g,  UPotential/N: %8g\n",
			    i,UKinetic/NumberOfParticles,UPotential/NumberOfParticles);
	    break;
    } else {
	if(UPotential<Uold) 
		EMStep *= 1.5;
	else 
		EMStep *= 0.5;
         Uold=UPotential;
      }
  }
  
  //The function PositionNONPDB is used to calculate the MSD
  for (j=0; j<NumberOfParticles; j++) {
	  PositionsNONPDB[j].x=Positions[j].x;
	  PositionsNONPDB[j].y=Positions[j].y;
	  PositionsNONPDB[j].z=Positions[j].z;
  }
  
  
 } 
